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ABSTRACT 


We consider three attributes of an individual that are critical in determining the temporal dynamics of 
pandemic influenza: social activity, proneness to infection, and proneness to shed virus and spread infec- 
tion. These attributes differ by individual, resulting in a heterogeneous population. We develop discrete- 
time models that depict the evolution of the disease in the presence of such population heterogeneity. For 
every individual, the value for each of the three describing attributes is viewed as an experimental value 
of a continuous random variable. The methodology is simple yet general, extending more traditional dis- 
crete compartmental models that depict population heterogeneity. Illustrative numerical examples show 
how individuals who have much larger-than-average values for one or more of the attributes drive the 
influenza wave, especially in the early generations of the pandemic. This heterogeneity-driven pandemic 
physics carries important policy implications. We conclude by using contact data in four European coun- 


tries to demonstrate empirical uses of our model. 


© 2012 Elsevier B.V. All rights reserved. 


1. Introduction 


Pandemic influenza is an ongoing threat to communities large 
and small. Despite a large body of research dealing with mathemat- 
ical modeling of influenza spread, there is still no widely accepted 
way to inform decision making in the case of a global event such 
as the 2009-2010 H1N1 “swine flu” pandemic. The H1N1 influenza 
outbreak demonstrated that limited information and short time 
allowances during a pandemic make it difficult to model accurately 
the spread of influenza and to make appropriate and effective deci- 
sions with regard to public policy and vaccine distribution. 

As a pandemic evolves, decision makers receive aggregate sta- 
tistics in the form of the number of people reporting to physicians 
with flu-like symptoms, number of related hospital admissions, 
number of flu-related deaths, and number of vaccinations adminis- 
tered. Yet, aggregate statistics hide the fact that early transmission 
and propagation of the disease are driven largely by particular seg- 
ments of the population: (1) those who are highly active in daily 
face-to-face encounters; (2) those who are overly prone to become 
infected given exposure; and/or (3) those who shed virus and 
spread the disease much more than average. Any person can be 
characterized along a spectrum of these three attributes: social 
activity, proneness to infection, and proneness to shed virus and 
spread infection. Those who are at the ‘right-hand-tails’ of one or 
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more of these distributional attributes play a significant role in 
the early spread of the disease. Such individuals, due to early infec- 
tion, drop out of the susceptible population near the middle and al- 
most certainly by the end of the outbreak. 

The correlation between these types of heterogeneity also plays 
a role in the speed of infection spread aside from the usual trans- 
missibility parameters like Ro. It is reasonable to assume that peo- 
ple who are most susceptible to infection are also those that are 
most likely to spread it to others; the most socially active people 
combine these two attributes. However, a negative correlation be- 
tween susceptibility and infectiousness will cause the infection to 
spread much slower and be sustained for a longer period of time 
within a community. To understand the dynamics of flu spread, 
or the spread of any human-to-human infectious disease, one can- 
not ignore such population heterogeneities. 

In this paper, we generalize the heterogeneous-population, 
influenza-spread model of Larson [15]. We improve analytical trac- 
tability by incorporating a more general, continuous form of popu- 
lation heterogeneity. The generalized model allows for a more 
realistic heterogeneous “class” selection, which in turn increases 
the accuracy of the model and yet simplifies its representation. 


2. Relevant modeling literature 


There is an extensive literature on modeling pandemic influ- 
enza spread. The simplest and most common models for influenza 
spread follow some variant of the S-I-R compartmental approach 
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where each person in the population is susceptible, infected, recov- 
ered, or dead. Such models are often used in a homogeneous set- 
ting, where all people in a single compartment behave identically 
and there is random mixing [1,14]. This approach may be consid- 
ered incomplete because it ignores population heterogeneity, the 
visible form of which is manifested as “super-spreaders”. The ex- 
treme cases of super-spreaders such as Typhoid Mary in the case 
of typhoid fever in the US [29], or Liu Janlun during the SARS out- 
break [18,26], bring with them the possibility of one particularly 
active person sparking an outbreak even if the rest of the popula- 
tion is relatively weak at transmitting the infection [16,22]. More 
typical are the populations where some members of the population 
are inherently more active at spreading the disease than others. 
Literature on heterogeneity consistently shows that this inherent 
heterogeneity causes a significant impact on the dynamics of the 
outbreak and even whether such an outbreak occurs at all 
[1,15,21]. Nigmatulina and Larson showed that people who have 
a large number of daily contacts are the “drivers” of the infection. 
These people tend to get infected early in the outbreak and drive 
the propagation of the flu in the early stages [15,21]. 

Heterogeneous compartmental models are also developed in 
the literature. The general idea is to divide the population into a fi- 
nite (usually small) number of discrete departments. People in 
each compartment are assumed to behave identically, with full in- 
tra-compartmental random mixing [13,23,30]. This approach in- 
volves a large number of differential equations with fixed 
Markovian assumptions, and needs to be solved numerically to 
compute the model-predicted average trends of the epidemic. In 
discrete-time heterogeneous models, a next-generation matrix is 
a useful tool in tracking an epidemic [9]. The next-generation ma- 
trix takes into account the anticipated inter-group contact rates, 
and uses them to calculate the number of infected individuals in 
each heterogeneous class of the population as the outbreak pro- 
gresses. It is used in modeling early-stage epidemic growth and 
in estimating flu parameters, such as the basic reproductive number, 
Ro, defined as ‘‘the expected number of secondary cases produced 
in a completely susceptible population by a ‘typical’ infected indi- 
vidual during her entire period of infectiousness” [10]. 

Employing an approach quite different from the traditional 
compartmental models, in this paper we expand on the heteroge- 
neous modeling by Larson [15]. Larson introduced simple discrete- 
time modeling that took into account population heterogeneity 
without the need for the data that forms the next-generation ma- 
trix. Instead, he used the idea of proportional mixing, where each 
person is likely to get infected in proportion to his or her contact 
rate, as opposed to a posited specific rate of infection between each 
pair of classes [15,23]. 

In this paper we generalize Larson’s approach to eliminate the 
need for a finite number of discrete classes of statistically identical 
individuals. Instead, we introduce a continuous distribution for all 
the parameters in question, in essence employing an infinite 
number of classes. Eventually we deal with all three attributes 
introduced above: social activity, proneness to infection, and 
proneness to spread infection. Initially we focus on contact rates, 
the measure of social activity. The model relies on just a few equa- 
tions that define the state of infection at a given time. We present 
numerical examples showing how various segments of the popula- 
tion are affected throughout the outbreak and the importance of 
determining a good starting distribution for contact rates. We also 
discuss the calculation of Ro, and particularly the effects of hetero- 
geneity on R,»- the average number of secondary infections caused 
by a typical infectious individual at time t from the start of the 
outbreak. We conclude by using contact data from four European 
countries to demonstrate the empirical uses of our model. In addi- 
tion to generalizing the traditional compartmental model by using 
a spectrum of heterogeneity as opposed to arbitrary classes, our 


model introduces closed-form equations to describe the state of 
the epidemic in a community at different points in time. These 
equations can be used to mathematically prove results regarding 
the changing characteristics of the population and the progression 
of the virus within the community. We provide some such proofs 
of in Section 3.3. 


3. Generalized model 


Larson employs a discrete-time model where the unit of time is 
a generation of influenza, or a “day”. A generation for the purposes 
of this model is the period of infection during which a person is 
infectious and actively interacting in society. Larson’s model as- 
sumes a finite number of classes of people, differing by their activ- 
ity level. Specifically, 


4c = Poisson rate of human contacts per day of an individual in 
class C. 

p = probability that a susceptible person becomes infected after 
a contact with an infectious individual. 


During each day i, a person may be classified as susceptible, in- 
fected, or immune. A person’s state may change on day i+ 1, where 
each susceptible person independently may become infected or re- 
main susceptible, an infected person recovers and becomes im- 
mune, and each immune person remains immune. For simplicity 
here, we assume that there are no influenza-related deaths during 
this outbreak. 

We generalize on this model by making the Poisson contact rate 
2a random variable, and assuming a distribution of 2 over the en- 
tire population. For a given person in the population, we sample 
from the j-distribution and find 4= 9. Then, that person has a 
number of potentially infectious daily contacts n, selected from a 
Poisson distribution with mean Jo. Here 49 corresponds to the indi- 
vidual’s “class”. We introduce the following notation: 


fd) population distribution of /, the individual Poisson rate of 
contacts per day! 
fA) distribution of , the individual Poisson rate of contacts per 


day i, in the infectious population on day. That is, if we 
were to select a single person, at random, from the infec- 
tious population,f/(4) would represent the distribution 
corresponding to that person’s daily contact rate, A 


f(A) distribution of 1, the individual Poisson rate of contacts per 
day, in the susceptible population on day i. That is, if we 
were to select a single person at random from the suscep- 
tible population, f* (A) would represent the density func- 
tion corresponding to that person’s daily contact rate, 4.7 

Ni number of infectious individuals on day i 

NS number of susceptible individuals on day i 

N total number of people in a population 


Once we know the initial conditions, we can model a typical 
epidemic curve and approximate the social-contact characteristics 
of the population by sequentially applying steps similar to those 
used by Larson: 


1 We assume here that f(4) is a continuous probability density function. The 
arguments in this section may also be applied to discrete distributions, for which f(/) 
represents the probability that the Poisson contact rate equals 7. We will use f(/) and 
the word “distribution” in both contexts to refer to either the probability density 
function or the probability mass function and use summations instead of integrals 
where appropriate. 

? Some individuals may change their contact rates during their infectious 
generation, e.g. by staying at home while symptomatic. We discuss the possibility 
of such behavioral changes in Section 3.4.3. 
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Day i: 

We start with the known values of all our quantities of interest 
during day i: N}, Ni, f§(4) and f!(J). 

We now show how to find all the necessary values of 
Niwa Nisa Sia and fi, (A). 

As per Larson, let 

fh; = the day i probability that a random person’s next interaction 
is with an infected person. 

We now consider the total expected number of contacts that all 
infectious people make on day i: 


| ” Nlafi(aaa = NI | * pl(ayda = NIEL(A). 
0 0 


The mean total number of day i contacts by infectious people is the 
product of the average number of contacts by an infectious person 
and the number of infectious people in the population. 

Likewise, the mean total number of human contacts each day is: 


[ ” Naf(ada = N | ” af (ada = NE(2). 
0 0 


and, with no behavioral changes or deaths, is assumed to be con- 
stant throughout the outbreak. 

Therefore, the probability that any interaction is with an infec- 
tious person, is simply the ratio: 


NiE;(A) 
Bi= "NEG: (1) 


We next define the day-specific infection probability: 
pi(I/2) = probability that a susceptible individual with contact 
rate 2 becomes infected on day i which Larson showed to be: 


pi(l/2) =1— eh, (2) 


We now have all the necessary tools to find N},,,Niv1; 3 4(4) and 
te 1(A): 

NG =f PUIA)Ne fe (Ada = ni fo P(IA)fP (Ada (3) 
Niwa = Ni — Niys- (4) 


The 4-distributions over different segments of the population are 
slightly more complex. The distribution of contact rates of infec- 
tious people changes on day i+ 1 from what it was on day i. On 
day i+1, it is the distribution of contact rates of those people 
who were susceptible and became infected on day i: 

1 ,(4) =f(Alindividual was infected on day i) 

= PAA) (5) 

© Jo PHAR) 


Now for the susceptibles: their day i+ 1 contact-rate distribution is 
that of those who were susceptible on day and who were not in- 
fected on dayit+1: 


= f(A|individual was not infected on day i) 


_ (1 — p,(I|4)) (a) 
Jo 1 = pill) ARAL 


These formulae describe entirely the state of the epidemic in the 
community on day i. By starting with some known boundary values 
of parameters and iteratively applying these rules over time, we ob- 
tain a model of complete epidemic dynamics within a continuous 
heterogeneous context. 


Aa (A) 


(6) 


Model summary 


Starting with the following quantities on day i: N?, Nj,f>(A) 
and f(A) 


Compute the intermediate quantities: 


pi(l|a) = 1 — eA 


alae 4 the relevant quantities on day: i+ 1 
= =N? Io Dj (| A)fp (Add 


5 S T 
ne, = Nj - ce 
PAL: 
A 
aA )F paaasoai D(I|A) ae 
(1—p,(Il4) i 


S a) ET PIMA A) 
ia (4) Jo Q-pidlangs (aya 


3.1. Numerical examples 


3.1.1. Uniform distribution 

We first model a population of a small hypothetical community 
of people. As a boundary condition, consider fo(4) to be uniform be- 
tween 0 and 40. Iteratively running our model, results are summa- 
rized in Figs. 1 and 2 and Table 1. Additional details are provided in 
Appendix C. 

Talking through the example, on day 0, we introduced infec- 
tious individuals into the population, each having a 2 value drawn 
from the uniform distribution over 0-40. However, on day 1, the 
distribution of contact rates of infectious persons resembles a 
straight line with slope 3. That is, if on day 1 we were to select a 
person at random from the infectious population, that person 
would be much more likely to have a high daily contact rate, than 
a low one. In fact, the expected contact rate by an infectious person 
on day 1 is 26.7 (Table 1). Note, also, that it is impossible for a per- 
son from the infectious population to have a daily contact rate of 0 
because a person with no contacts has no chance of being infected. 
The infected population becomes “high-activity” heavy. This is 
consistent with Larson’s model and supports the conclusion that 
high-activity people are the drivers of the infection at the begin- 
ning of the outbreak. As the outbreak continues, however, the 
high-activity tail dips down, as high-activity individuals become 
infected early, become immune and can no longer be infected. Note 
that, on the contrary, the distribution of contact rates in the sus- 
ceptible population becomes more and more “Jow-activity” heavy. 
On the last day of the outbreak, the average contact rate in the sus- 
ceptible population is 10.4, almost half of the original average con- 
tact rate of 20. 


3.1.2. Poisson distribution 

We now consider the same community with different initial 
conditions. We propose a more realistic, Poisson distribution for 
f(4) with mean 20. Firstly, we are now dealing with a discrete dis- 
tribution for 2 as opposed to a continuous density function. Fig. 3 
describes the progression of the outbreak with some notable differ- 
ences from the first example. 

Here too, we see that the characteristics of the subsections of 
the population change as the outbreak continues, although the ef- 
fect is much smaller. On the first day, the expected contact rate in 
the infectious population is 21, a slight increase from 20(Table 2). 

Toward the end, the expected contact rates in different subsec- 
tions of the population change slightly, but generally the distribu- 
tions remain similar to the original contact rate distribution, 
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Epidemic Curve R(t) 
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Fig. 1. Under an initial uniform distribution of daily contact rates, the epidemic curve and R(t) of the outbreak. R(t) is computed by dividing the number of infectious people at 


time t+ 1 by the number of infectious people at time t. 


Contact Rate PDFs for Infecteds 
0.1 1 r ; 


Contact Rate PDFs for Susceptibles 
0.1 : + : 


Fig. 2. Contact rate characteristics for the infectious and susceptible subsections of the population at different times during the outbreak. 


Table 1 
Outbreak statistics. The initial distribution of 4 is uniform between 0 and 40. 


Outbreak statistics 


Duration of outbreak 17 days 
Total number of infections 7168 
EI (A) 26.7 
Fi7(2) 18.2 
Et7() 10.4 


(Fig. 4). This behavior is noticeably different from the changing 
distributions in the uniform distribution case. 


3.2. Sensitivity to original distribution 


The two numerical examples show that the original distribution 
of contact rates plays a major role in the evolution of all four of our 


relevant quantities. The total expected number of infections in the 
outbreak is 10% larger in the case of a Poisson original distribution 
than in the case of a uniform distribution. In the case of the uni- 
form distribution, the susceptible and the infectious subsections 
of the populations have extremely different contact rate character- 
istics. The susceptible population is low-activity heavy, while the 
contact rate distribution for the infectious population is a concave 
function with low probability at the left tail and decreasing proba- 
bility at the right. 


3.2.1. Calculation of Ro 

The differences between the two cases are particularly evident 
in the calculation of the basic reproductive number, Ro. Recall that 
Ro is defined as the expected number of secondary infections 
caused by a typical infected individual in a fully susceptible popu- 
lation [10]. On day O, since almost everybody is susceptible, we can 
assume that we are in a fully susceptible population. There is some 
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Fig. 3. The epidemic curve and R(t) of the outbreak. The initial distribution is Poisson with mean 20. 


Table 2 
Outbreak statistics. The initial distribution of 4 is Poisson with mean 20. 


Outbreak statistics 


Duration of outbreak 21 days 
Total number of infections 7865 
E\ (A) 21.0 
E4o(A) 19.5 
ES (A) 18.4 


uncertainty as to the exact manner of calculating Ro. The cause of 
ambiguity is the term “typical”. Suppose we pick a person at ran- 
dom from our infectious population on day O. In both the case of 
the uniform and the Poisson original distributions, the average 
contact rate for a random person, call her patient zero, would be 


Contact Rate PMFs for Infecteds 


0.1 r r CT CT 
Day 0 
0.09 + Day 1/7 
Day 10 
0.08 + Day 20/1 


probability 
oO 
fo} 
a 


20. If the Poisson contact rate is 20, the average number of patient 
zero’s contacts is also 20. Since each contact has .1 probability of 
spreading infection, this interpretation of the definition would give 
us an Rg of 2. If Ro is the same for both original distributions of con- 
tact rates, we would expect the total number of infections to be 
approximately the same, but they are not. 

Instead, we recommend the logic followed by network models 
in the established epidemiological literature [9,10,19], where Ro 
is calculated as an expected value of a distribution that is weighted 
by the magnitude of the contact rate. That is, a person in the infec- 
tious population with a high contact rate should be weighed pro- 
portionally higher, as that person is more likely to acquire and 
then spread the infection than a person with a low contact rate. 
In other words, our “typical” person would no longer be selected 
at random, but according to a distribution fr(2). Such a distribution 
will no longer be sampling from the population of people, but rather 


Contact Rate PMFs for Susceptibles 
0.1 : v ; ; 


Day 0 
Day 1 |] 
Day 10 
Day 20]] 


probability 
oO 
lo} 
oa 


Fig. 4. Contact rate characteristics for the infectious and susceptible subsections of the population at different times during the outbreak. The initial distribution is Poisson 


with mean 20. 
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population of contacts. Consider each person dropping a piece of 
paper that is marked with that person’s 2 value on the floor during 
each and every contact. We then pick Rp by drawing a random 
piece of paper. This new distribution would favor high activity 
Spd 
EA (4) 

Note that this is exactly the formula used to describe the con- 
cept of time-average duration of an interval in a renewal process 
[12]. Using this scenario, we can calculate Ro. For the uniform dis- 
tribution case: 


40 We A402 
p dh = 
> | AOE) 30% 20 


people proportionally to their activity level, then f(A) = 


Ro 27. 


For the Poisson case: 


oo 92 
Py AOD 
Po= > Paaale 2007 21 


Note that even though the uniform distribution has a larger value of 
Ro, it actually results in fewer total infections! This phenomenon is 
caused by the fact that the uniform distribution has a higher vari- 
ance than a Poisson distribution with the same mean. In the uni- 
form case, there are fewer very high-activity people who are 
depleted much more quickly than in the Poisson case. Once the 
high-activity people are gone, the outbreak quickly subsides result- 
ing in a shorter total outbreak time and fewer total infections. It is 
important to note, also, that the effect of depletion of high-activity 
individuals from the population comes from the fact that high- 
activity people are more susceptible to infection; we show in Sec- 
tion 3.4.2, that the fact that these people are also more infectious 
does not contribute to this effect. Andersson and Britton [2] con- 
firm, using compartmental model methods, that epidemics are par- 
ticularly devastating when susceptibility does not vary, i.e. there is 
no early depletion. 

This is another reminder that an influenza model is a simplifica- 
tion of a real life problem that is too complex to describe fully. To 
avoid critical mistakes, healthcare decision-makers must be extre- 
mely cautious in choosing initial values of the model’s parameters 
such as the initial distribution of contact rates to avoid critical mis- 
takes. Many published heterogeneous modeling published papers 
describe deriving Ro in a manner consistent with the second ap- 
proach above, arguing that Rp is the mean of a distribution that 
is skewed toward high-activity individuals [9,10,19]. However, 
one needs to be careful when transferring this idea into simulation 
modeling. Often, we find that flu simulations are seeded by select- 
ing a few typical “patient zeros” from the population to enter the 
susceptible population in an infectious state. These typical people 
are selected to be representative of the overall general population, 
and their average contact rate is used to calibrate the simulation to 
known Rp parameters. “Patient zeros” that are selected in this way, 
result in an underestimate of the true value of Ro. 


3.3. Common properties 


Despite the fact that the starting distribution plays a significant 
role in the way an epidemic curve for the outbreak is calculated 
using this model, certain properties remain true regardless of ori- 
ginal distribution. 


Theorem 1 (Tilting of contact rate distributions in susceptible and 


infected populations). Starting with any initial conditions, f}(A) and 
f§ (A): for any (i> 0): 


fig=Gle @ -e FRM), (7) 


and 
FU=C le = EU), (8) 


where Ci and C are normalizing constants. Equations (7) and (8) im- 
ply that as the outbreak progresses, contact rate distributions become 
tilted distributions with tails decaying exponentially faster than in 
the original contact rate distributions. Note also, that the initial distri- 
bution of infectious individuals, f}(/), plays no role in the distributions 
for any generation i > 0. 


Proof. * 


We next introduce a result that will be useful in describing the 
behavior of distributions as the outbreak progresses in more detail. 


Lemma 1. For any constants 0<a<1, 0<b<1, and j4>0 the 
function: 


av = (GO) 


is monotonically decreasing and convex in 2. 


Theorem 2. Stochastic ordering of susceptible and infectious popula- 
tions. Let A} be the random variable distributed as f5(2), and similarly 
A} as f5(4). Then, as long as f(A), fi(2) are non-deterministic, for any 
i> 0: 


Ai,i<srA; for 2 > E?(A), (9) 
Af<spAi for 2 > E}(A), (10) 
Ai,,<srA? for 2 > E}(A), (11) 


where X<ST implies that X is strictly smaller than Y in the usual 
stochastic order [24]. Stochastic ordering of the various distributions 
as time goes on provides us with an analytical expression of our intu- 
ition. That is, according to Eq. (9), the contact rates in the susceptible 
population become “smaller” with time in the stochastic ordering sense. 
In other words, the susceptible population becomes more and more 
densely packed with inactive individuals as the high activity people be- 
come infected. Similarly with the infectious population. 


Theorem 3. Depletion of high-activity individuals in the susceptible 
and infectious populations. As long as f5(2),fi(4) are non-determinis- 
tic, for any i> 0: 


fo >fh for > EA), (12) 
fo<srfi, for 2 > E}(A), (13) 
fi>stfi., for 2 > Ei(A). (14) 
These are true regardless of initial contact rate distribution. Eq. (12) 


implies that after the first generation, the susceptible population be- 
comes increasingly low-activity heavy. In Eq. (14) we show that the 
distribution of infected people also becomes more low-activity heavy 
with time since high-activity people become depleted from the popula- 
tion. For illustration, see Fig. 5. We note that since we model infection 
spread in discrete time increments, the high-activity people get infected 
in large numbers within each generation. However, in a real-life 
scenario, these generations overlap and we do not need to “wait for a 
generation to end” for new people to get infected. As a result, the deple- 
tion of high activity individuals should occur even more rapidly in a 
continuous-time, real-world outbreak. Finally, we note that expectation 


3 All proofs in this section are presented in Appendix A. 
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Lemma 3 Illustration 


Eq9 Eq 10 


Susceptibles Day 10 
Susceptible Day 11 | 0.09 } 


l 
Infected Day 10 
Infected Day 11 


Susceptibles Day 10 
Infected Day 11 


probability 


0 10 20 30 40 
nr 


Fig. 5. Illustration to Theorem 3: the probability density functions of contact rates for the susceptible and infected populations on days 10 and 11. The initial distribution of 


contact rates was uniform, ranging from 0 to 40. 


of contact rate distributions in both the susceptible and infected subsec- 
tions of the population is decreasing with time. 


Corollary 1. Overall decrease in expected contact rates. For i> 0 and 
non-deterministic f§(A),f§(A): 


Ei (A) < Ei,4(4), (15) 
EA) <E4(). (16) 


3.4, Further generalizations 


The methodology developed above can be extended to include 
heterogeneity due to vulnerabilities of different groups as well as 
the effects of changing social behaviors during an ongoing 
pandemic. 


3.4.1. Mortality 

So far in our model, we have assumed that all infected people 
recover and rejoin the population within one generation. How- 
ever, any devastating pandemic influenza that we fear is likely 
to be deadly; some infected individuals are likely to die of the 
infection. Moreover, the likelihood of death from an infection will 
almost certainly not be uniform throughout the population. The 
Centers for Disease Control and Prevention (CDC) suggests that 
pregnant women and the elderly are at most risk from death dur- 
ing a seasonal influenza outbreak [6,17]. Meanwhile, one of the 
characteristics of a particularly dangerous novel influenza strain 
pandemic, is the heightened mortality rates of the people in their 
20s [27]. 

To add this layer of complexity, we add a probability of death 
after an infection, q. That is, an infectious person dies at the end 
of her infectious period with probability q, and rejoins the popula- 
tion with probability (1 — q). Moreover, this probability cannot be 
constant, as certain people are at more risk than others. We must 
impose a probability distribution on q that is not independent from 
Ae 

The quantities of interest become: 
fiZ.q) the joint distribution of 4, the Poisson rate of contacts per 

day and, q the probability of death given infection, in the 
entire population alive at the start of day i 


the joint distribution of 4, the Poisson rate of contacts per 
day, and q, the probability of death given infection, in the 
infectious population on day i 

the joint distribution of 4, the Poisson rate of contacts per 
day, and q, the probability of death given infection, in the 
susceptible population on day i 

Ni the number of people alive in the population at the start of 
day i 

as before 


fq 


mC) 


Ni,N; 


Replicating Eqs. (1)-(6): 
The total number of contacts made by infectious people on day i 
remains as before: 


[ [mn 


However, the total number of interactions made by the entire pop- 
ulation on day i, now needs to account for the mortality rate of the 


population: 
q)didq = N; [ [a aU 


[ [ Natiea 


We can now solve for f;: 


NiE;(A) 
hie NiEi(A)’ 


q)didq = Nj - - NiAfl(A, q)dadq = NjEi(A). 


q)didq = NiE;(2). 


(Ja) 


The dynamics of the spread of infection on day i remain unchanged 
as the deaths occur at the end of each generation; Eqs. (3)-(6) re- 
main unchanged except for the new distribution functions that ac- 
count for death probabilities. 


Nhe [ [ paanife@.aididg =i a [ paare.qaiad, 


(3a) 

Nia = Nj - Ne (4a) 
PULA), @) 

ted (A, q) lags (I) Naee qjda (Sa) 
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(1 — pA) APA.) 

1 (As | — L 6a 
val to (1 — p,(I|2)) 5(4,q)dA X62) 

We are left to find the effects of the deaths on day i. Consider the 

probability that a person with contact rate 2 and mortality probabil- 

ity q died in period i: 


P;(death|1,q) = P;(infection|A, q)q, 
NI 
AGN 


AOE 
Fi(A,q) — fi(A,q)Ni’ 


P;(infection|A, q) 


for all values of 2, q for which fi(1,q) > 0. 
So the probability that a person with these parameters survived 
day i, is then: 


1 AON 
fil, qQ)Ni 


We can then derive fi+1(4,q): 


‘ fi.qNt 
faa(.9) fi(A, q)Pi(no death|,, q) fU.g(1 — fates) (17) 
we ed P(no death) 1 — NiEa) 
Nj 


and finally, assuming that deaths are independent of each other: 


Nia =Ni- wif f° afi(2 


3.4.2. Variability in susceptibility and infectivity 

The second additional layer of complexity is also acknowledged 
by Larson. We have so far assumed that if a susceptible and infec- 
tious individual come into contact, regardless of their respective 
contact rates, the probability that a susceptible person becomes in- 
fected from the contact is p. However, particularly in recent times, 
there has been added concern about people who are inherently 
more susceptible to disease than others. For example, the recent 
H1N1 pandemic of 2009 was shown to be less likely to affect peo- 
ple over 65 [5]. The “Asian Flu” of 1957 reported similar tendencies 
most likely due to acquired immunity by the older generation. 
Similarly, some portion of the population includes “super-shed- 
ders”; these people are much more likely to spread the disease 
than the average person. These differences in people’s ability to 
shed and be infected by the flu virus, cause variability in the 
parameter p. 

Larson brings up this possibility by splitting p into two different 
parameters. We will change notation here to call these parameters 
rand s. 

r =the infectiousness index of an individual. r can take on any 
value between 0 and 1; it describes the level of infectiousness of 
an individual, with 0 implying a person does not shed any virus, 
and 1 implying a person is maximally infectious. 

s =the level of susceptibility of an individual. That is, the prob- 
ability that a susceptible individual gets infected after an interac- 
tion with an infectious individual with infectiousness index r= 1. 

Using this notation, if an infectious person with infectiousness 
index r and a susceptible person with susceptibility parameter s 
make contact on a given day, the probability that the susceptible 
person becomes infected is p=rx«s 

We alter our epidemic parameters to include r and s: 

fi(4,q,17,S) = joint distribution of 4,q,r, and s in the entire popula- 
tion alive at the start of day i. 


q)didq = Ni — NjEi(q). (18) 


f!(,q,1,8) = joint distribution of 2,q,r, and s in the infectious 
population on day i. 
f$(4,q",5) = joint distribution of 2,q,r, and s in the susceptible 


population on day i. 


N; = number of people alive in the population at the start of day 


Ni, NP as before. 

Once again we calculate equations analogous to Eqs. (1)-(6), 
taking extra care with notation: 

Br) = the probability that a random person’s next interaction is 
with an infected person having infectivity parameter r, on day i. 
Jo Jo So" NiAft( 


4,q,1,5)dadqds NiE.(Alr)fi(r) 


‘(r : 1b 
Ba) fo So So So’ NiAfi, 4, 1, s)dadqdrds NE(A) on 
Let: 

p,(I|2,q,S) = Probability that a susceptible person with 


susceptibility parameter s and contact-rate 
parameter 1, becomes infected on day i 
(note that the person’s value of r is irrelevant), (2b) 
Pan tees 
pilia,q,s) =1—e* Jo * 


See Larson [15] supplementary materials for derivation. 
Then we calculate the usual quantities: 


al 1 il oe) 
Nia = Ni f i | [ D(A, 4, SFP (4, g, 1, s)dAdqdrds, (3b) 


NL, = NS- Ni 


i4+1 i+1? 


pi(la,q, s)fS(A, 4, 7,5) 
(12,4, S)FS(A, q, 1, s)dadqdrds * 


is (A, q; r, Ss) 
Cie 


(1 = pill, 4.5) f?(4,4,1,8) 


fi4(4,4,1,8) 

mee fo fo i Jo A — pUl\4, 4,5) FP A.4GT ,s)dddqdrds’ 
(6b) 

TA, NI 
f(A,4,8)(1— Reread 
fir (2,9,1,8) ( _- ! (14b) 
1-Sp@ 
and 
Nia = Ni — NiEi(Q). (15b) 


3.4.3. Behavioral changes 

We have so far been working with a static model, where people 
behave in the same way every day and do not change their routine 
as influenza enters their community. However, previous outbreaks 
such as the 2002 SARS epidemic show evidence that people change 
their behavior throughout their outbreak even if their region has 
not yet been affected by the epidemic [21,25]. These behavioral 
changes have been shown to make an impact on the spread of 
infection. A good model, therefore, should take into account the 
fact that people will try to limit their daily activity in the event 
of a pandemic scare [21]. 

This can be easily incorporated our model by introducing a 
transformation function G1) that represents the effect of behav- 
ioral changes in 4 on day i. If G is invertible, the distribution func- 
tions can be easily determined at the end of each day. 

On day i consider the quantities f(/,q,r,s),f! (A,q,1,s) and 
f? (4,q,1,s) to be the relevant distributions computed without tak- 
ing into account behavioral interventions. Then: 


fi(4,4,1,8) =i (Gil4),4.1,8), 
f1240,8) =f (Gl4).4.0,5), 
FAG.1,8) =F (G4), 4,7,5). 
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For example, we can assume that each day, all individuals reduce 
their number of contacts by a factor of 4 . Then G,(A) = 44, and 
the distributions are changed accordingly. This is a simple example 
of how the model can be used to reflect dynamic behavioral changes 
during the outbreak. The model can be extended further and made 
arbitrary complex by using non-deterministic or state-based func- 
tions for partial compliance with government programs. 


3.5. European data 


We conclude with an example based on research by Mossong 
et al. In this section, we will use the word ‘day’ to refer to an actual 
real-time day, and the word ‘generation’ to refer to a generation of 
influenza. The Mossong group conducted a thorough study of con- 
tacts by participants in eight European countries. Participants were 
asked to record their daily contacts, “defined as either skin-to-skin 
contact such as a kiss or handshake... or a two-way conversation 
with three or more words in the physical presence of another per- 
son” [20]. The information from the participants’ diaries was 
weighed to match the demographics of participating countries. 
The group published distributions of daily contacts by individuals 
from each country. We ran our model on the populations of four 
of those countries, namely, Belgium, Great Britain, Germany, and 
Poland, to demonstrate empirical use of our model. 

To incorporate this information into our model, we note that 
these data were selected by sampling distinct persons and noting 
the number of contacts that person had in one day. This is not ex- 
actly f() as we have defined it for two reasons. Firstly, f(/) is the 
distribution of contact rates within a generation, not one day. To 
correct for this, we assumed that a generation of the flu is approx- 
imately 2.5 days [31,28], and extrapolated the available data 
accordingly. Secondly, f(,) is the distribution of contact rates within 


heterogeneity based solely on contact rates, and none of the other 
possible heterogeneous parameters (see Tables 3 and 4). 

Table 5 summarizes the results of our model for each of the four 
countries. Figs. 6 and 7 further illustrate the epidemic curves and 
the characteristics of the population before and after the 
outbreaks. 

It is not surprising with all transmissibility parameters being 
equal, that the Ro as well as the total number of infections caused 
by the outbreak is correlated with the mean contact rate in the 
population. It is also worth noting, that as suggested in our model 
analysis, the variance in the contact rates in the country is corre- 
lated with change in the susceptible population characteristics that 
occurred throughout the outbreak, see Table 6. Recall, also, that 
variance plays a role in the Ro calculations, and as a result while 
Belgium and Great Britain have similar average contact rates in 
their respective populations, the values of Ro for the two countries 
differ. 

Higher variance implies that a country has a combination of 
highly active people and relatively non-active people in its popula- 
tion. During the first few generations of the outbreak, the high 
activity super-spreaders become infected and immune, while the 
susceptible population becomes dominated by low-activity 
individuals. 


4. Discussion 


Current policy literature gives little consideration to changing 
characteristics of the population seen throughout the outbreak. 
However, an argument can be made, that heterogeneity should 
play a significant role in certain policy decisions. Consider for 


a population, and not the distribution of the number of daily con- Table 5 
tacts by an individual, which is the quantity being sampled. Recall The results of running the model using initial parameters described in Tables 3 and 4. 
that the number of daily contacts by a person is Poisson distributed Belgium Great Britain Germany Poland 
with mean Each sample point in the empirical data is a value Saale 
picked from this Poisson distribution. However, we note that this Ry 1.98 1.66 133 254 
empirical distribution is the maximum-likelihood estimator for Total 4,620,000 25,780,000 15,788,000 23,787,000 
f(A) in a specific country and we used it to estimate f(,) for our infections (44%) (41%) (19%) (62%) 
model. If the data from Mossong constitute a representative sam- Duration 103 days 118 days 258 days 83 days 
ple of the population, then the daily contact distribution should 
be equivalent to f(). We used the same basic parameters for each 
country, so that the only differences between countries affecting 
the epidemic curve were their respective sizes and their contact 
rate distributions. Finally, note that in this example we model ModeiGensraled Epidemia Gunes 
15 T : T T T 
Belgium 
Great Britain 
Table 3 Germany 
Initial global parameters used in model instance. Poland 
A mo} 
Basic parameters a) 
Oo 
qi) .002 (estimated H1N1 mortality rate) [7] £ 
P 04 c 
Ni 10 2 
oO 
SD 
Qa 
{o) 
a 
xe) 
xe 
Table 4 
Initial country dependent parameters used in model. 
Belgium Great Germany Poland 
Britain 
Country-dependent parameters 
Population 10,423,493 62,348,447 82,282,988 38,463,689 0) 50 100 150 200 250 300 
[8] [8] [8] [8] Day 
Average daily 11.84 11.74 7.95 16.31 


contact rate 


Fig. 6. Epidemic curves for the 4 countries as percentage of the population infected 
on a given day. 
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Contact Rate Distributions in Susceptible Population 


Belgium 
0.25 
HE Prior to Outbreak 
0.2 HE After Outbreak 
2 0.15 
= 
[sj 
ne} 
2g 0.1 
[om 
0.05 , | | | 
| PRR RRE EE 
0 10 20 30 40 50 
Contacts 
Germany 
0.25 
HE Prior to Outbreak 
0.2 HE After Outbreak 
2 0.15 
= 
[sj 
Ke} 
g 0.1 
Qa 
0.05 | | 
ge . i | a 
0 10 20 30 40 50 
Contacts 


Great Britain 


0.2 1 
Prior to Outbreak 
| 
= After Outbreak 
0.15 
0.1 4 
0.05 ] | | | 
TREE 
0 10 20 30 40 50 
Contacts 
Poland 
0.2 
Prior to Outbreak 
=a 
0.15 i] After Outbreak 
0.1 | 
- [ | | L 
0 L i L 1 1 L | 
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Contacts 


Fig. 7. The contact rate distributions in the susceptible population of each country before the beginning and after the end of the outbreak. 


Table 6 
The variance of a region’s contact rate distribution compared to the total change in its 
population characteristics. 


Belgium — Great Germany Poland 
Britain 
Contact rate variance and population characteristic change 
Contact rate variance 9.85 7.67 6.26 11.45 
Mean square population 6.6e—005 4.8e—005 1.0e—005 6.6e—004 


change 


example, vaccine administration. In a population with little heter- 
ogeneity and low contact rate variance, the characteristics do not 
change significantly throughout the outbreak, and vaccines 
administered early in an outbreak reach similar classes of people 
as vaccines administered later on. In a high-variance population, 
vaccines administered early will be more likely given to high-activ- 
ity people, whereas vaccines administered later on will be particu- 
larly ineffective, as most of the recipients will already be low- 
activity individuals. We must not forget, however, that other con- 
siderations must play a role in policy decisions. Researchers have 
suggested that it is often preferable to vaccinate high-risk individ- 
uals like pregnant women and people with comorbidities [17]. 
However, mathematical models such as this one may provide in- 
sights into when and why activity-focused vaccination policies 
may be preferable. 

The model we provide in this paper is generalizes existing het- 
erogeneous compartmental models. The continuous spectrum of 
heterogeneous properties allows modelers a natural alternative 
to discrete compartments with identical individuals and propor- 
tional inter-compartmental mixing. Inclusion of heterogeneity in 
population modeling is especially important in policy decisions 
such as occur in designing vaccine allocation plans and message- 


targeting in public awareness campaigns. The framework we pro- 
vide not only simplifies modeling heterogeneous populations, but 
also provides insight into changes of population attributes over 
time. The closed form equations describing evolving population 
characteristics should invite further mathematical analysis of 
epidemics. 
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Appendix A. Proofs 


Theorem 1. Starting with any initial conditions, f}(A) and f(A); for 
any i> 0: 
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flay=cfe "6 7B" lay, (A7) 
pac] e "| F(a, (A8) 


Where Cj and C; are normalizing constants. 


Proof. We proceed by induction. Suppose that for some i > 0, 


-p > 
fla=Cle i 18a) 
Then 
i-1 
NANPA) 1 Beye ee aa Poe 
| Z Di( Cc 1 e ApB; e j=0 ‘Ss yi 
1 ( ) Jo P (ASS (A)dA int ) Ko ) 
; > h nS ; 
=G.,[e © -e  [fp(r). 


Since is a constant [> p,(I|A)fP(A). Similarly: 


(= pulayfia) cs ce-iony fg 25" | 
mHO= Fa paapaa Me™ le |e) 
SF 
=Gale ™ |fe(2) 


Now for the base case; note that since there were no infected people 
prior to day 0, B_; =0: 
Po(llA)fo(A) 1 
fi@az = Gl 
(= Te pydtas@d — “"" 
= Cy (1 — e 7%) (e-P) f(A) 


e-/PH)F5(2) 


Ce fo (2) = Cie?) 


Lemma 1. For any constants 0<a< 1,0<b< 1, the function: 


#0)= (ere) 


is monotonically decreasing and convex in 1 for 1 > 0. 


Proof. We first prove convexity. 
Let Ab = a and # = y and h(a) = (Ge eT) 7) .Then g(A) is convex in 2 
if and only if h(a) is convex in a. Differentiating: 


Ph e-*((e2 + e*)e% + (—€2% + 2e% — 1) y? + (2e* 
doz e3% — 3e24 4 3e% — 1 
The denominator e?% — 3e?% + 3e% — 1 =(e* — 1)? > O since «> 0. 

Consider now the numerator. The numerator evaluated at «=0 
is 0. We want to show that it is positive for all other «, and 
therefore the numerator is positive for all « > 0. Let e* = w and let’s 
factor out the positive term e~*’. We are left to examine: 


2e2*)y e2% e*) 


2 


were BY — (y+ 1)? + (2p? + 2y— Iu. 
Differentiating with respect to pL: 


dk ; 7 
Tied yer? + (1+ py? — 2(y + 1)? + (2y? + 2y - 1). 
Evaluating at = 1: 
dk 
A (2+ y)+1+y-2(y+1) + (2 + 2y- 1) =0. 
LU u=1 
Differentiating again: 
dK ; 
ae (1+ (2+ yw + (1+ yyw! - 2(y 41) 
dK 
qe, =~ A+N2+y)+A+ 7-27 +1" =0, 
u=1 
dK 


qe VE ee 
Since all quantities in the above equation are positive, the third 
derivative x of is positive for all 4. > 1. Note also, « that and its first 
and second derivatives are 0 at w= 1, which implies that the first 
and second derivatives of « as well = x itself must be positive for 
all >1.. This in turn implies that £4 is positive for all «>0O and 
the function g is convex for all 1 >0. 

To prove the monotonic descent property, we note that h(«) is 
monotonically decreasing if and only if g(4) is monotonically 
decreasing: 


on 


dh e-*"(e*/+* + (1 — e*)y — e*) 


da (e* — 1)? 
erie =1) filer) ) 
(e*—1)? <0 


For all « > 0. Consequently,g() is monotically decreasing. 


Theorem 2. Stochastic ordering of susceptible and_ infectious 
populations. 

Let A? be the random variable distributed as f°, and similarly A? as 
FPA). 


Then, as long as f§(A),fi(4) are non-deterministic, for any i> 0: 


A’,,<srA? for 4 > E;(A), (A.9) 
Ai<srAj,, for 2 > E(A), (A.10) 
Ai<st4; for 2 > Ei(A), (A.11) 


where X < 5; implies that X is strictly smaller than Y in the usual sto- 
chastic order. 


Proof. Take any i> 0: 
(= PAA _ p53) 
x1 — pi(llA))f7(4) 


(e~/P6i) 
E(e-?h) 


£02) en) 


= (2) 
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Since f£;>0,e-?% is strictly monotonically decreasing in /. 


Therefore: 
aA) (e7/Phi) 
B(4) — Ei(e-?h) 


is strictly monotonically decreasing in 2 implying that A? i;1 iS smal- 
ler than 4? in likelihood ratio order, which in turn implies Eq. (9). 
See Shaked and Shanthikumar for further detail [24]: 


pi(llAyastf?(4) _ ps,,. (1—e-**) 
i41 (A) ior) \FS f (A A) S 
Do Pi(llA)P (A) E3(1 — e-Phi) | 
Here (1 — e~’?) is strictly monotonically increasing, so: 
RA _ Bd — eA) 
1) (leh) 


Is strictly monotonically decreasing in 2, implying Eq. (10). 
Finally, 


i-2 i-1 
DS ee DOr) ft  -s 
flay=cfe Pe 18" | (2) = A (em — 1912), 


i 


ie) 
asi (4) 
ae) = sy’ 


Furthermore: 


is 
7 ahi) (1 — ePhi) 
(EPA — 1)" HEi(! (1—e Phi a) 


e’PPi-1 —1) 


(1 — e- Phi) 
E?(1 — e~Phi) 


(1-e“PAi) 
— py) hee) 
! B (4 (1- yo 
e’PBi-1 1) 
From the definition of £; , and the fact that i> 0; £j_; > 0 and f;>0. 


Ap bij 
From Lemma 1, we know that (He is monotonically 
decreasing, so: 


na (A) =F (A) 


(1-e7/Phi) 


ha (A) (Ge 
fia B(; (1—e7 PB; +) 


e“Phi-1 —1) 


is monotonically decreasing, and Eq. (11) is proved. 


FA)>R,U) forrz> Bd, (A.12) 
FP) < fia) for 2 > E?(), (A.13) 
fi) >fis@) fora > EA), (A.14) 


regardless of initial contact rate distribution. 


Proof. Take any i> 0. 
By Jensen’s inequality, since e~’?4: is a strictly convex function, 
eH MPBi < E5(e~/Phi) 
Furthermore, since e~/PAi 
> EF(A). 
e7 PB <e 
eae 2. ; 
—Ap 
(4) = SPA De “hy L < f5(4), proving Eq. (12). 


Now for Eq. (13), note that since f§(4) is non-deterministic and 
4 > 0,E}(4) > 0. We are interested in 4 > E$(4) > 0. 

Since (1—e-’PF:) is strictly concave for 2>0 , we again use 
Jensen’s inequality, to show that: 


is strictly decreasing, whenever 


E})Pbi < E}(e7Phi) , and using the derivation from 


1—e Fire > E5(1 —e 
E5(4)1 — e- 


(1 — e-/PAi) 
E}(1 — e-*Phi) 


7h) 


whenever 2 > PB > 1 — e-FAPBi S E5(1 — e-/Phi): 


.1(4) RFA > fP(A). 


Finally, we use Jensen’s inequality for the third time for the last 
equation. Recall from Lemma 1, ey that is monotonically 
decreasing and convex for 4 > 0. 

Whenever 4 > E!(2): 


(1 — e-Pbi) (1- e-Fi()Pbi) ((1 — e-Phi) 
(a —1) S (e8()Pbi-1 1) < £; (ePhia — 1))? 


and 


iad) <fild)- 


Corollary 1. For i> 0, and f§(2),f}(4) non deterministic: 


F(A) < Ei, (A), 
Ei(A) < Ei, (A). 


cate 
_ 


Proof. The proof follows immediately from the stochastic ordering 
of A? and Aj. 


Appendix B. Supplementary material 


Supplementary data associated with this article can be found, in 
the online version, at doi:10.1016/j.ejor.2012.01.027. 
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